*** Remotely incorrect application replication ***
*** Date: 10 August 2022 *******
****************************************
***Naive models: No Misclassification***
****************************************

**Spec 1: Full Sample, Ejido FEs**

clear all
set more off
cap log cl

*cd "C:\Users\03638881\Dropbox\Satellite data measurement\"
*glo data "C:\Users\03638881\Dropbox\Satellite data measurement\Data\"
*glo outputs "C:\Users\03638881\Dropbox\Satellite data measurement\Latex\tables\"
*adopath + "C:\Users\03638881\Dropbox\Satellite data measurement\Do Files\"

local 	sel_user		2		// 1 - Dan, 2 - Jen
	
	if 		`sel_user' 		== 1 { 
		global 	user		"C:\Users\03638881\Dropbox\Satellite data measurement\"
	}	
	
	if 		`sel_user' 		== 2 {
		global user 		"C:\Users\alixgarj\Dropbox\Satellite data measurement\"
	}
	

	
	global	data			"${user}Data\"
	global 	outputs			"${user}Latex\tables\"
	global 	graphs	 		"${user}Latex\graphs"
	global 	dofiles			"${user}Do files\"

adopath + "${user}Do Files\"


log using "${outputs}MeasurementErrori2iDataJAG1", replace
use "${data}PESi2iPanelDataShortAugust2021.dta", clear
qui {
	loc spec = 1
	keep if inejido==1
	loc fe "EjidoCode"
	loc cluster "EjidoCode"
	sort APEJNCNID EjidoCode yr
	xtset `fe'

	*define vars to be used
	loc y "binary_defor"

	xi i.yr
	loc x0 "_Iyr_2002 _Iyr_2003 _Iyr_2004 _Iyr_2005 _Iyr_2006 _Iyr_2007 _Iyr_2008 _Iyr_2009 _Iyr_2010 _Iyr_2011 _Iyr_2012 _Iyr_2013 _Iyr_2014"
	loc x1 "beneL1 everbene elev_met_mean slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind"
	loc x2 "beneL1 everbene c.TPApAreaHa##c.path_row_counts_2000_2015_L7 c.slope_deg_mean##c.path_row_counts_2000_2015_L7 dist_any_road_mt dist_5000_city_km elev_met_mean pctTC2000 per_maj_ind"

	*rescale vars*
	replace TPApAreaHa=TPApAreaHa/100
	replace dist_any_road_mt=dist_any_road_mt/1000
	replace dist_5000_city_km=dist_5000_city_km/100
	replace slope_deg_mean=slope_deg_mean/100
	replace elev_met_mean=elev_met_mean/100

	g z1=path_row_counts_2000_2015_L7
	g z2=path_row_counts_2000_2015_L7*slope_deg_mean 
	g z3=path_row_counts_2000_2015_L7*TPApAreaHa
	lab var path_row_counts_2000_2015_L7 "Cloud-free scenes, L7"
	lab var z2 "Cloud-free scenes, L7 x Avg Slope"
	lab var z3 "Cloud-free scenes, L7 x Area"
	lab var dist_5000_city_km "Distance to city with $>$ 5,000 people"

	loc z1 "path_row_counts_2000_2015_L7"
	loc z2 "z2 z3"

	*create means for CRE models*
	loc mx1 " "
	loc mz1 " "
	loc mz2 " "
	foreach x in `x1' {
		bys `fe': egen m`x'=mean(`x')
		loc mx1 = "`mx1' m`x'"
	}
	foreach x in `z1' {
		bys `fe': egen m`x'=mean(`x')
		loc mz1 = "`mz1' m`x'"
	}
	foreach x in `z2' {
		bys `fe': egen m`x'=mean(`x')
		loc mz2 = "`mz2' m`x'"
	}
}

loc i=1
xtreg `y' `x0' `x1', fe i(`fe') vce(cl `cluster')
eststo col`i'
loc i=`i'+1
xtreg `y' `x0' `x1' `z1' `z2' , fe i(`fe') vce(cl `cluster')
eststo col`i'
loc i=`i'+1
logit `y' `x0' `x1' `mx1', vce(cl `cluster') 
loc ll = e(ll)
eststo col`i': margins, dydx(`x1') predict(pr) post gen(l_`spec'_)
eststo col`i', add(logL `ll')


loc i=`i'+1
logit `y' `x0' `x2' `mx1' `mz1' `mz2', vce(cl `cluster')
loc ll = e(ll)
eststo col`i': margins, dydx(`x1' `z1') predict(pr) post gen(ahl_`spec'_)
eststo col`i', add(logL `ll')

loc i=`i'+1

loc logL = -999999999

forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap scobit `y' `x0' `x1' `mx1', vce(cl `cluster') constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo col`i': margins, dydx(`x1') predict(pr) post gen(s_`h'_`spec'_)
				eststo col`i', add(alpha `a')
				eststo col`i', add(logL `ll')
			}
		}
		
	}	
	

loc i=`i'+1

loc logL = -999999999


forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap	scobit `y' `x0' `x2'  `mx1' `mz2', vce(cl `cluster') constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo col`i': margins, dydx(`x1' `z1') predict(pr) post gen(ahs_`h'_`spec'_)
				eststo col`i', add(alpha `a')
				eststo col`i', add(logL `ll')
			}
		}
		
	}	
	

loc i=`i'+1


esttab col*  ///
	using "${outputs}i2iresults1.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" "" "" "") ///
	 keep(`x1' `z1' `z2') nonumbers  replace  fragment  ///
	 scalar(alpha logL)  prehead({ \begin{tabular}{lcccccc}    ///
	 \hline \hline ///
	 & 		  & Ad Hoc &  & Ad Hoc &  & Ad Hoc  \\  ///
	 & FE-LPM & FE-LPM & CRE Logit & CRE Logit & CRE Scobit & CRE Scobit\\  ///
	 & (1) & (2) & (3) & (4) & (5) & (6)  \\ ) ///
	 prefoot( \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate estimator. ///
	  Standard errors are clustered at the municipality level. Marginal effects evaluated at sample means are displayed for the logit and scobit models. ///
	  Fixed effects are at the Ejido level. Time fixed effects also included in all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

estimates clear 

keep `y' `x1' `z1' `z2'  yr muni_cve i2itpapid EjidoCode APEJNCNID ahl* ahs* l_* s_*
save "${data}PESi2iPanelData_resultsJAG.dta", replace


************************************************************************************************

**Spec 2: Full Sample, polygon FEs**

clear all
use "${data}PESi2iPanelDataShortAugust2021.dta", clear
qui {
	loc spec = 2
	keep if inejido==1
	loc fe "APEJNCNID"
	loc cluster "EjidoCode"
	sort APEJNCNID EjidoCode yr
	xtset `fe'

	*define vars to be used
	loc y "binary_defor"

	xi i.yr
	loc x0 "_Iyr_2002 _Iyr_2003 _Iyr_2004 _Iyr_2005 _Iyr_2006 _Iyr_2007 _Iyr_2008 _Iyr_2009 _Iyr_2010 _Iyr_2011 _Iyr_2012 _Iyr_2013 _Iyr_2014"
	loc x1 "beneL1"
	loc x2 "beneL1 c.TPApAreaHa##c.path_row_counts_2000_2015_L7 c.slope_deg_mean##c.path_row_counts_2000_2015_L7"

	g z1=path_row_counts_2000_2015_L7
	g z2=path_row_counts_2000_2015_L7*slope_deg_mean 
	g z3=path_row_counts_2000_2015_L7*TPApAreaHa
	lab var path_row_counts_2000_2015_L7 "Cloud-free scenes, L7"
	lab var z2 "Cloud-free scenes, L7 x Avg Slope"
	lab var z3 "Cloud-free scenes, L7 x Area"

	loc z1 "path_row_counts_2000_2015_L7"
	loc z2 "z2 z3"

	*create means for CRE models*
	loc mx1 " "
	loc mz1 " "
	loc mz2 " "
	foreach x in `x1' {
		bys `fe': egen m`x'=mean(`x')
		loc mx1 = "`mx1' m`x'"
	}
	foreach x in `z1' {
		bys `fe': egen m`x'=mean(`x')
		loc mz1 = "`mz1' m`x'"
	}
	foreach x in `z2' {
		bys `fe': egen m`x'=mean(`x')
		loc mz2 = "`mz2' m`x'"
	}
}

loc i=1
xtreg `y' `x0' `x1', fe i(`fe') vce(cl `cluster')
eststo col`i'
loc i=`i'+1
xtreg `y' `x0' `x1' `z1' `z2' , fe i(`fe') vce(cl `cluster')
eststo col`i'
loc i=`i'+1
logit `y' `x0' `x1' `mx1', vce(cl `cluster') 
loc ll = e(ll)
eststo col`i': margins, dydx(`x1') predict(pr) post gen(l_`spec'_)
eststo col`i', add(logL `ll')

loc i=`i'+1
logit `y' `x0' `x2' `mx1' `mz1' `mz2', vce(cl `cluster')
loc ll = e(ll)
eststo col`i': margins, dydx(`x1' `z1') predict(pr) post gen(ahl_`spec'_)
eststo col`i', add(logL `ll')

loc i=`i'+1

loc logL = -999999999


forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap	scobit `y' `x0' `x1'  `mx1', vce(cl `cluster') constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				loc ll = e(ll)
				eststo col`i': margins, dydx(`x1') predict(pr) post gen(s_`h'_`spec'_)
				eststo col`i', add(alpha `a')
				eststo col`i', add(logL `ll')

			}
		}
		
	}	
	

loc i=`i'+1

loc logL = -999999999


forval a=  0.2(0.05)0.95 {
    loc h=`a'*100
	loc aa = ln(`a')
	constraint 1 /lnalpha = `aa'
	cap	scobit `y' `x0' `x2'  `mx1' `mz2', vce(cl `cluster') constraint(1) 
		if e(converged) == 1 {
			if e(ll) > `logL' {
				cap drop xb
				predict xb, xb	
				loc logL = e(ll)
				eststo col`i': margins, dydx(`x1' `z1') predict(pr) post gen(ahs_`h'_`spec'_)
				eststo col`i', add(alpha `a')
				eststo col`i', add(logL `ll')

			}
		}
		
	}	
	

loc i=`i'+1


esttab col*  ///
	using "${outputs}i2iresults2.tex", ///
	se nonotes  style(tex)  b(%12.3f) se(%12.3f)  noobs ///
	starlevels(* 0.10 ** 0.05 *** 0.01) label mlabels("" "" "" "" "" ""  ) ///
	 keep(`x1' `z1' `z2') nonumbers  replace  fragment  ///
	  scalar(alpha logL) prehead({ \begin{tabular}{lcccccc}    ///
	 \hline \hline ///
	 & 		  & Ad Hoc &  & Ad Hoc & & Ad Hoc \\  ///
	 & FE-LPM & FE-LPM & CRE Logit & CRE Logit & CRE Scobit & CRE Scobit \\  ///
	 & (1) & (2) & (3) & (4) & (5) & (6) \\ ) ///
	 prefoot( \\ ) ///
	  postfoot(\hline \hline \end{tabular} } \begin{tablenotes}[para,flushleft] \footnotesize{Column headers indicate estimator. ///
	  Standard errors are clustered at the municipality level. Marginal effects evaluated at sample means are displayed for the logit and scobit models. ///
	  Fixed effects are at the polygon level. Time fixed effects also included in all models. * p $<$.10, ** p$<$ .05, *** p$<$.01.} \end{tablenotes} )

estimates clear 


keep yr EjidoCode APEJNCNID ahl* ahs* l_* s_*
merge 1:1 yr EjidoCode APEJNCNID using "${data}PESi2iPanelData_resultsJAG.dta"
tab _m
drop _m
save "${data}PESi2iPanelData_resultsJAG.dta", replace


use "${data}PESi2iPanelData_resultsJAG.dta", clear

#delimit ;
order yr APEJNCNID EjidoCode i2itpapid TPApAreaHa per_maj_ind dist_any_road_mt dist_5000_city_km muni_cve elev_met_mean slope_deg_mean binary_defor beneL1 
pctTC2000 everbene path_row_counts_2000_2015_L7 z2 z3 l_1* ahl_1* s_20_1* ahs_20_1* s_20_2* 
ahs_20_2* ;
#delimit cr

*Logit, Scobit
loc i=1
foreach v in beneL1 everbene elev_met_mean slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind {
	foreach s in 1 {
		ren l_`s'_`i' l_`s'_`v'
		foreach a in 20 {
			ren s_`a'_`s'_`i' s_`a'_`s'_`v'
		}
	}
	loc i=`i'+1
}
loc i=1
foreach v in beneL1 {
	foreach s in 2 {
		ren l_`s'_`i' l_`s'_`v'
		foreach a in 20 {
			ren s_`a'_`s'_`i' s_`a'_`s'_`v'
		}
	}
	loc i=`i'+1
}

*Ad Hoc Logit, Scobit
loc i=1
foreach v in beneL1 everbene elev_met_mean path_row_counts_2000_2015_L7 slope_deg_mean dist_any_road_mt dist_5000_city_km TPApAreaHa pctTC2000 per_maj_ind {
	loc vv = "`v'"
	if ("`v'"=="path_row_counts_2000_2015_L7") loc vv = "path"
	foreach s in 1 {
		ren ahl_`s'_`i' ahl_`s'_`vv'
		foreach a in 20 {
			ren ahs_`a'_`s'_`i' ahs_`a'_`s'_`vv'
		}
	}
	loc i=`i'+1
}
loc i=1
foreach v in beneL1 path_row_counts_2000_2015_L7 {
	loc vv = "`v'"
	if ("`v'"=="path_row_counts_2000_2015_L7") loc vv = "path"	
	foreach s in 2 {
		ren ahl_`s'_`i' ahl_`s'_`vv'
		foreach a in 20 {
			ren ahs_`a'_`s'_`i' ahs_`a'_`s'_`vv'
		}
	}
	loc i=`i'+1
}

save "${data}PESi2iPanelData_results_cleanJAG.dta", replace






log close